#include "monte-carlo.h"

/* Retorna um número aleatório entre 0 e RAND_MAX */
double mc_random_number(drand48_data *status)
{
    double val;
    
    drand48_r(status, &val);
    
    return val;
}

/* Calcula dois números aleatórios e verifica se eles estão dentro de 1/4 da circunferência */
int mc_is_inside_circle(drand48_data *status)
{
    double x = mc_random_number(status);
    double y = mc_random_number(status);
    
    return (x*x + y*y <= 1);
}

/* Calcula a quantidade de pontos dentro da circunferência de ITS pontos aleatórios */
unsigned int mc_number_of_points_in_circle(drand48_data *status)
{
    unsigned int i, j;
    
    for(i = j = 0; i < POINTS; i++)
        if(mc_is_inside_circle(status)) j++;
    return j;
}

/* Calcula o valor de pi */
void mc_pi_value(mpf_t pi, mpz_t n, mpz_t m)
{
    mpf_t aux;
    mpf_init(aux);

    mpf_set_z(aux, n);
    /* Calcula o valor de Pi usando a fórmula pi = 4*m/n */
    mpf_set_z(pi, m); mpf_mul_ui(pi, pi, 4); mpf_div(pi, pi, aux);

    mpf_clear(aux);
}

/* Inicializa os parâmetros de Monte Carlo */
void mc_init_param(mc_param *p)
{
    int i;

    mpf_inits(p->pi[0], p->pi[1], NULL);
    mpz_init(p->n);
    p->stop = B_FALSE;

    for(i = 0; i < NTHREADS; i++)
    {
        srand48_r(time(NULL), &(p->c[i].status)); /* Semeia o algoritmo para calcular números aleatórios */
        p->c[i].i = &(p->i); /* Aponta para o contador da estrutura "pai" */
        mpz_inits(p->c[i].m[0], p->c[i].m[1], NULL);
    }
}

/* Limpa os parâmetros de Monte Carlo */
void mc_clear_param(mc_param *p)
{
    int i;

    mpf_clears(p->pi[0], p->pi[1], NULL);
    mpz_clear(p->n);

    for(i = 0; i < NTHREADS; i++)
        mpz_clears(p->c[i].m[0], p->c[i].m[1], NULL);
}

/* Algoritmo sequencial de Monte Carlo */
void monte_carlo(mpf_t pi, long *it)
{
    mpf_t old_pi; /* Variável para armazenar os valores antigos de Pi */
    mpz_t m, n;
    int i;
    drand48_data status;

    srand48_r(time(NULL), &status); /* Semeia o algoritmo para calcular números aleatórios */

    mpf_init(old_pi);
    mpz_inits(m, n, NULL);

    *it = 0;
    for(i = 0; i < MAX_ITS; i++)
    {
//        printf("%d\n", i);

        /* Incrementa a quantidade total de pontos */
        mpz_add_ui(n, n, POINTS);
        /* Incrementa a quantidade de pontos dentro da circunferência */
        mpz_add_ui(m, m, mc_number_of_points_in_circle(&status));

        /* Seta a variável old_pi com o valor de pi */
        mpf_set(old_pi, pi);

        /* Calcula o valor de pi */
        mc_pi_value(pi, n, m);

        /* Compara old_pi com pi e, se houver convergência, pare */
        if(has_converged(old_pi, pi)) break;
    }

    *it = i; /* Número de iterações utilizadas para calcular pontos aleatórios */

    mpz_clears(m, n, NULL);
    mpf_clear(old_pi);
}

/* Algoritmo concorrente de Monte Carlo */
void monte_carlo_conc(mpf_t pi, long *it)
{
    pthread_t thread[NTHREADS+1]; /* Serão NTHREADS para calcular valores de m e 1 thread para calcular Pi */
    mc_param p;
    int i, j;

    mc_init_param(&p);

    for(i = 0; i < MAX_ITS_CONC && !p.stop; i++)
    {
        p.i = i;
//        printf("%d\n", i);

        /* Criação das NTHREADS threads que calcularão valores de m */
        for(j = 0; j < NTHREADS; j++)
            pthread_create(&thread[j], NULL, mct_calculate_points, &p.c[j]);

        /* Criação da thread que calculará o valor de Pi e verificará convergência */
        pthread_create(&thread[NTHREADS], NULL, mct_calculate_pi, &p);

        /* Espera todas as threads acabar para sincronizá-las */
        for(j = 0; j <= NTHREADS; j++)
            pthread_join(thread[j], NULL);
    }

    *it = i*NTHREADS; /* Número de iterações utilizadas para calcular pontos aleatórios */
    mpf_set(pi, p.pi[INDEX2(i, 1)]);
    mc_clear_param(&p);
}

/* Função da thread para calcular números aleatórios e fazer a contagem de m */
void *mct_calculate_points(void *arg)
{
    mc_count *c = (mc_count *) arg;
    int i0, i1;

    /* i0 terá o valor de 0 quando i for par e 1 quando i for ímpar */
    i0 = INDEX2(*(c->i), 0);
    i1 = INDEX2(*(c->i), 1);

    /* Incrementa o valor de m */
    mpz_add_ui(c->m[i1], c->m[i0], mc_number_of_points_in_circle(&(c->status)));

//    printf("m = ");
//    mpz_out_str(stdout, 10, c->m[i1]);
//    printf("\n");

    return NULL;
}

/* Função da thread que calcula o Pi e verifica convergência */
void *mct_calculate_pi(void *arg)
{
    mc_param *p = (mc_param *) arg;
    mpz_t m;
    int i, i0, i1;

    /* Não é feito o cálulo do Pi quando i for 0, pois não temos valores de m para tal */
    if(p->i == 0) return NULL;

    /* Atribuição dos índices i0 (índice do valor antigo) e i1 (índice do valor novo) */
    i0 = INDEX2(p->i-1, 0);
    i1 = INDEX2(p->i-1, 1);

    /* Quantidade de pontos totais utilizados é a quantidade de threads vezes ITS */
    mpz_add_ui(p->n, p->n, NTHREADS * POINTS);

    /* Realiza a somatória total dos pontos dentro da circunferência */
    mpz_init(m);
    for(i = 0; i < NTHREADS; i++)
        mpz_add(m, m, p->c[i].m[i1]);

    /* Calcula o valor de Pi */
    mc_pi_value(p->pi[i1], p->n, m);

    /* Se houve convergência, a variável stop é setada como true */
    if(has_converged(p->pi[i0], p->pi[i1])) p->stop = B_TRUE;

    mpz_clear(m);

    return NULL;
}
